#include <errno.h>
#include <gmp.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>

/* Default number of iterations */
#define ITERATIONS 25

/* Maximum number of iterations */
#define MAX_ITERATIONS 100

/* Default minimum precision, in bits */
#define PRECISION 35000000

static void *thread_x(void *arg);
static void *thread_y(void *arg);
static void *thread_p(void *arg);

/* Variables shared between multiple threads */
static mpf_t x[MAX_ITERATIONS+1],
        sqrt_x[MAX_ITERATIONS+1],
        inv_sqrt_x[MAX_ITERATIONS+1],
        y[MAX_ITERATIONS+1],
        p0, p;

/* One mutex per shared variable per iteration */
static pthread_mutex_t mut_x[MAX_ITERATIONS+1],
        mut_sqrt_x[MAX_ITERATIONS+1],
        mut_inv_sqrt_x[MAX_ITERATIONS+1],
        mut_y[MAX_ITERATIONS+1];

/*
 * This thread calculates all values of x, storing them into x[].
 *
 * At first, we lock all mut_x[], mut_sqrt_x[] and mut_inv_sqrt_x[] mutexes and
 * then unlock them as new values are calculated so that other threads can use
 * them.
 */
void *
thread_x(void *arg)
{
    int i, iterations;

    iterations = *((int *) arg);
    for (i = 0; i <= iterations; ++i) {
        if (i > 0) {
            pthread_mutex_lock(&mut_x[i]);
            pthread_mutex_lock(&mut_sqrt_x[i]);
        }
        pthread_mutex_lock(&mut_inv_sqrt_x[i]);
    }

    /* fprintf(stderr, "thread_x: calculating x[0]\n"); */
    mpf_init_set_d(x[0], 2.0);
    mpf_sqrt(x[0], x[0]);
    pthread_mutex_unlock(&mut_x[0]);

    for (i = 1; i <= iterations; ++i) {
        /* fprintf(stderr, "thread_x: calculating x[%d]\n", i); */

        mpf_init(sqrt_x[i-1]);
        mpf_sqrt(sqrt_x[i-1], x[i-1]);
        pthread_mutex_unlock(&mut_sqrt_x[i-1]);

        mpf_init_set(inv_sqrt_x[i-1], sqrt_x[i-1]);
        mpf_ui_div(inv_sqrt_x[i-1], 1, sqrt_x[i-1]);
        pthread_mutex_unlock(&mut_inv_sqrt_x[i-1]);

        mpf_init(x[i]);
        mpf_add(x[i], sqrt_x[i-1], inv_sqrt_x[i-1]);
        mpf_div_ui(x[i], x[i], 2);
        pthread_mutex_unlock(&mut_x[i]);
    }

    return NULL;
}

/*
 * This thread calculates all values of y, storing them into y[]. Each loop
 * waits until the next x is available for use.
 *
 * At first, we lock all mut_y[] mutexes and then unlock them as new values are
 * calculated so that other threads can use them. Whenever a variable
 * calculated by another thread is needed, we first trylock its mutex because
 * if the mutex is unlocked (retval != EBUSY), that variable is ready to be
 * used.
 */
void *
thread_y(void *arg)
{
    int i, iterations;
    mpf_t tmp_num, tmp_den;

    iterations = *((int *) arg);
    for (i = 1; i <= iterations; ++i)
        pthread_mutex_lock(&mut_y[i]);
    mpf_inits(tmp_num, tmp_den, NULL);

    /* fprintf(stderr, "thread_y: calculating y[1]\n"); */
    if (pthread_mutex_trylock(&mut_sqrt_x[0]) == EBUSY)
        pthread_mutex_lock(&mut_sqrt_x[0]);
    mpf_init_set(y[1], sqrt_x[0]);
    pthread_mutex_unlock(&mut_y[1]);

    for (i = 2; i <= iterations; ++i) {
        /* fprintf(stderr, "thread_y: calculating y[%d]\n", i); */

        if (pthread_mutex_trylock(&mut_sqrt_x[i-1]) == EBUSY)
            pthread_mutex_lock(&mut_sqrt_x[i-1]);
        mpf_mul(tmp_num, y[i-1], sqrt_x[i-1]);

        if (pthread_mutex_trylock(&mut_inv_sqrt_x[i-1]) == EBUSY)
            pthread_mutex_lock(&mut_inv_sqrt_x[i-1]);
        mpf_add(tmp_num, tmp_num, inv_sqrt_x[i-1]);

        mpf_add_ui(tmp_den, y[i-1], 1);

        mpf_init(y[i]);
        mpf_div(y[i], tmp_num, tmp_den);
        pthread_mutex_unlock(&mut_y[i]);
    }

    mpf_clears(tmp_num, tmp_den, NULL);

    return NULL;
}

/*
 * This thread calculates all values of p, but only stores the current and
 * previous calculated values. Whenever a variable calculated by another thread
 * is needed, we first trylock its mutex because  if the mutex is unlocked
 * (retval != EBUSY), that variable is ready to be used.
 */
void *
thread_p(void *arg)
{
    int i, iterations;
    mpf_t tmp_num, tmp_den;

    iterations = *((int *) arg);
    mpf_inits(p0, p, tmp_num, tmp_den, NULL);

    /* fprintf(stderr, "thread_p: calculating p[0]\n"); */
    mpf_set_ui(p0, 2);
    if (pthread_mutex_trylock(&mut_x[0]) == EBUSY)
        pthread_mutex_lock(&mut_x[0]);
    mpf_add(p0, p0, x[0]);

    for (i = 1; i <= iterations; ++i) {
        /* fprintf(stderr, "thread_p: calculating p[%d]\n", i); */

        if (pthread_mutex_trylock(&mut_x[i]) == EBUSY)
            pthread_mutex_lock(&mut_x[i]);
        mpf_add_ui(tmp_num, x[i], 1);

        if (pthread_mutex_trylock(&mut_y[i]) == EBUSY)
            pthread_mutex_lock(&mut_y[i]);
        mpf_add_ui(tmp_den, y[i], 1);

        mpf_div(tmp_num, tmp_num, tmp_den);
        mpf_mul(p, p0, tmp_num);

        /* Store p into p0 for the next iteration */
        mpf_set(p0, p);
    }

    mpf_clears(tmp_num, tmp_den, NULL);

    return NULL;
}

int
main(int argc, char *argv[])
{
    int i, iterations, precision;
    char *p_str;
    mp_exp_t exp;
    pthread_t thread_x_handle, thread_y_handle, thread_p_handle;

    iterations = ITERATIONS;
    /* The number of iterations; if given but is not valid or too big, abort */
    if (argc > 1 && (sscanf(argv[1], "%d", &iterations) != 1 || iterations > MAX_ITERATIONS))
        exit(1);

    precision = PRECISION;
    /* The precision; if given but is not valid, abort */
    if (argc > 2 && scanf(argv[2], "%d", &precision) != 1)
        exit(1);

    mpf_set_default_prec(precision);

    for (i = 0; i <= iterations; ++i) {
        pthread_mutex_init(&mut_x[i], NULL);
        pthread_mutex_init(&mut_sqrt_x[i], NULL);
        pthread_mutex_init(&mut_inv_sqrt_x[i], NULL);
        pthread_mutex_init(&mut_y[i], NULL);
    }

    /*
     * These are locked before all threads begin because most of them will need
     * these values for initialization, and therefore must wait for them to be
     * calculated
     */
    pthread_mutex_lock(&mut_x[0]);
    pthread_mutex_lock(&mut_sqrt_x[0]);

    pthread_create(&thread_x_handle, NULL, thread_x, &iterations);
    pthread_create(&thread_y_handle, NULL, thread_y, &iterations);
    pthread_create(&thread_p_handle, NULL, thread_p, &iterations);

    /*
     * Wait for this thread as it will end only when we have our Pi
     * approximation
     */
    pthread_join(thread_p_handle, NULL);

    p_str = mpf_get_str(NULL, &exp, 10, 0, p);
    printf("%s\n", p_str + exp);

    /* Cleanup */
    free(p_str);
    for (i = 0; i <= iterations; ++i) {
        pthread_mutex_destroy(&mut_x[i]);
        pthread_mutex_destroy(&mut_sqrt_x[i]);
        pthread_mutex_destroy(&mut_inv_sqrt_x[i]);
        pthread_mutex_destroy(&mut_y[i]);

        mpf_clears(x[i], y[i], sqrt_x[i], inv_sqrt_x[i], NULL);
    }
    mpf_clears(p0, p, NULL);

    return 0;
}
